Clustering

Hierarchical clustering et Kmeans

Anne Badel, Frédéric Guyon & Jacques van Helden

2020-03-11

Les données

Les données dans l’ordinateur (1)

Les iris de Fisher

Ces données sont un classique des méthodes d’apprentissage Fisher

Les données dans l’ordinateur (2)

  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4

Les données dans l’ordinateur (2)

  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4

! : convention différente en RNA-seq

Représentons ces données : une fleur (1)

mes.iris[1,]
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2

Comment représenter cette fleur ?

Dans quel espace de réprésentation ?

Représentons ces données : une fleur (2)

plot(mes.iris[1,1:2])

Dans le plan, un point de coordonnées :

représenté par un vecteur \(v2 = (\) 5.1 \(,\) 3.5\()\) dans \(\mathbb{R}^2\)

Représentons ces données : une fleur (3)

Dans l’espace, un point de coordonnées :

représenté par un vecteur \(v3 = (\) 5.1 \(,\) 3.5 \(,\) 1.4\()\) dans \(\mathbb{R}^3\)

Représentons ces données : toutes les fleurs (4)

= un nuage de points dans un espace à 4 dimensions

= PAS de représentation possible (pour l’instant)

Représentons ces données : une variable à la fois (1)

Représentons ces données : deux variables à la fois (2)

Il faut tenir compte de toutes les dimensions

c’est à dire de toutes les variables à notre disposition

Clustering et classification (termes anglais)

On a une information sur nos données

Clustering : on cherche à mettre en évidence des groupes dans les données

Clustering et classification (termes anglais)

On a une information sur nos données

Clustering : on cherche à mettre en évidence des groupes dans les données

Classification :

Clustering

données simulées : y a-t-il des groupes ?

données simulées : y a-t-il des groupes ?

Géométrie et distances (1)

On considère les données comme des points de \(\mathbb{R}^n\)

\(\mathbb{R}^n\) : espace Euclidien à \(n\) dimensions, où

Géométrie et distances (2)

On considère les données comme des points de \(R^n\) (*)

Géométrie et distances (3)

Sur la base d’une distance (souvent euclidienne)

Distances

Définition d’une distance : fonction positive de deux variables

  1. \(d(x,y) \ge 0\)
  2. \(d(x,y) = d(y,x)\)
  3. \(d(x,y) = 0 \Longleftrightarrow x = y\)
  4. Inégalité triangulaire : \(d(x,z) \le\) d(x,y)+d(y,z)

Si 1,2,3 : dissimilarité

Distances utilisées dans R (1)

Distances utilisées dans R (2)

Autres distances non géométriques (pour information)

Utilisées en bio-informatique:

\[d("BONJOUR", "BONSOIR")=2\]

Distances plus classiques en génomique

Il existe d’autres mesures de distances, plus ou moins adaptées à chaque problématique :

Ne sont pas des distances, mais indices de dissimilarité :

Note : lors du TP, sur les données d’expression RNA-seq, nous utiliserons le coefficient de corrélation de Spearman et la distance dérivée, \(d_c = 1-r\)

Avec R (1) : distance entre deux individus

1.04 2.34 1.59 4.12 3.69
2.13 2.68 3.18 2.84 3.68

distance euclidienne : 4.98

distance de manhattan = 13.48

Avec R (2) : distance entre individus d’un nuage de points

mat.iris <- mes.iris[sample(1:150, 5),]
print(dist(mat.iris), digits = 2)
      5   33   43   72
33 0.56               
43 0.73 1.22          
72 3.13 3.19 3.40     
49 0.33 0.42 1.05 2.98
cor.mat.iris <- cor(t(mat.iris))
print(as.dist(1 - cor.mat.iris), digits = 2)
          5       33       43       72
33 0.002904                           
43 0.000078 0.002457                  
72 0.219978 0.252287 0.217709         
49 0.000318 0.004885 0.000415 0.205493

Avec R (3) : distance entre variables décrivant le nuage de points

cor.mat.iris <- cor(mat.iris)
print(as.dist(1 - cor.mat.iris), digits = 2)
             Sepal.Length Sepal.Width Petal.Length
Sepal.Width        1.3204                         
Petal.Length       0.1408      1.7215             
Petal.Width        0.1996      1.8069       0.0095

Distances entre groupes (1)

Distances entre groupes (2)

\[D(C_1,C_2) = \min_{i \in C_1, j \in C_2} D(x_i, x_j)\]

\[D(C_1,C_2) = \max_{i \in C_1, j \in C_2} D(x_i, x_j)\]

Distances entre groupes (3)

\[D(C_1,C_2) = \frac{1}{N_1 N_2} \sum_{i \in C_1, j \in C_2} D(x_i, x_j)\]

\(d^2(C_i,C_j) = I_{intra}(C_i \cup C_j)-I_{intra}(C_i)-I_{intra}(C_j)\)

\(D(C_1,C_2) = \sqrt{\frac{N_1N_2}{N_1 + N_2}} \| m_1 -m_2 \|\)

Distances entre groupes (4)

Les données

Ces données sont un classique des méthodes d’apprentissage

Dans un premier temps, regardons les données

dim(mes.iris)
[1] 150   4
head(mes.iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
str(mes.iris)
'data.frame':   150 obs. of  4 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
summary(mes.iris)
  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  

Visualisation des données

On peut ensuite essayer de visualiser les données

plot(mes.iris)

Préparation des données (1) : variables de variance nulle

iris.var <- apply(mes.iris, 2, var)
kable(iris.var, digits = 3, col.names = "Variance")
Variance
Sepal.Length 0.686
Sepal.Width 0.190
Petal.Length 3.116
Petal.Width 0.581
sum(apply(mes.iris, 2, var) == 0)
[1] 0

Préparation des données (2) : “Normalisation”

Afin de pouvoir considérer que toutes les variables sont à la même échelle, il est parfois nécessaire de standardiser les données.

mes.iris.centre <- scale(mes.iris, center = TRUE, scale = FALSE)
mes.iris.scaled <- scale(mes.iris, center = TRUE, scale = TRUE)

On peut visuellement regarder l’effet de la standardisation :

standardisation par la moyenne et la variance

par(mfrow = c(1,2))
par(mar = c(7, 4.1, 4.1, 1.1)) # adapt margin sizes for the labels
boxplot(mes.iris, main = "Raw data", las = 2)
boxplot(mes.iris.scaled, main = "scaled", las = 2)

par(mar = c(5.1, 4.1, 4.1, 2.1)) # Restore original margin sizes
par(mfrow = c(1,1))

standardisation par la médiane

iris.mediane <- apply(mes.iris, 2, median)
mes.iris.scaled.mediane <- sweep(mes.iris, 2, iris.mediane)
par(mfrow = c(1,2))
par(mar = c(7, 4.1, 4.1, 1.1)) # adapt margin sizes for the labels
boxplot(mes.iris.centre, main = "scaled, moyenne", las = 2)
boxplot(mes.iris.scaled.mediane, main = "scaled, médiane", las = 2)

par(mar = c(5.1, 4.1, 4.1, 2.1)) # Restore original margin sizes
par(mfrow = c(1,1))

par l’intervalle inter-quartile (IQR)

iris.iqr <- apply(mes.iris, 2, IQR)
mat.iris.iqr <- matrix(rep(iris.iqr, each = nrow(mes.iris)), ncol = 4)
mes.iris.scaled.iqr <- mes.iris / mat.iris.iqr
par(mfrow = c(1,2))
par(mar = c(7, 4.1, 4.1, 1.1)) # adapt margin sizes for the labels
boxplot(mes.iris, main = "raw data", las = 2)
boxplot(mes.iris.scaled.iqr, main = "scaled, IQR", las = 2)

par(mar = c(5.1, 4.1, 4.1, 2.1)) # Restore original margin sizes
par(mfrow = c(1,1))

La matrice de distances

Nous utilisons ici la distance euclidienne sur données normalisées.

La classification hiérarchique : principe

classification hiérarchique : mettre en évidence des liens hiérachiques entre les individus

Notion importante, cf distances

L’algorithme : étape 1

au départ

identification des individus les plus proches

construction du dendrogramme

étape j :

calcul des nouveaux représentants ‘BE’ et ‘CD’

calcul des distances de l’individu restant ‘A’ aux points moyens

A est plus proche de …

dendrogramme

pour finir

dendrogramme final

Je ne fais pas attention à ce que je fais …

… c’est à dire aux options des fonctions dist() et hclust()

par(mfrow = c(2, 1))
plot(iris.hclust, hang = -1, cex = 0.5, main = "Données brutes")
plot(iris.scale.hclust, hang = -1, cex = 0.5, main = "Normalisées")

En utilisant une autre métrique

En utilisant un autre critère d’aggrégation

En conclusion

Les heatmap

pheatmap::pheatmap(mes.iris)

pheatmap::pheatmap(mes.iris, scale = "column")

pheatmap::pheatmap(mes.iris, scale = "row")

Les k-means

Les individus dans le plan

=> faire apparaitres des classes / des clusters

L’algorithme

étape 1 :

Choix des centres provisoires

Calcul des distances aux centres provisoires

Affectation à un cluster

Calcul des nouveaux centres de classes

Etape j :

Fin :

Arrêt :

Un premier k-means en 5 groupes

iris.scale.kmeans5 <- kmeans(mes.iris.scaled, center=5)
iris.scale.kmeans5
K-means clustering with 5 clusters of sizes 53, 17, 47, 26, 7

Cluster means:
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1  -0.05005221 -0.88042696    0.3465767   0.2805873
2  -1.39493454 -0.05056417   -1.3357516  -1.3187694
3   1.13217737  0.08812645    0.9928284   1.0141287
4  -0.93018708  1.05972279   -1.2791040  -1.2202264
5  -0.38011687  2.26106915   -1.2952890  -1.1986013

Clustering vector:
  [1] 4 2 2 2 4 5 4 4 2 2 4 4 2 2 5 5 5 4 5 4 4 4 4 4 4 2 4 4 4 2 2 4 5 5 2 2 4 4 2 4 4 2 2 4 4 2 4 2 4 4 3 3 3 1 1 1 3 1 1 1 1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 3 3 3 1 1 1 1 1 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 3 3 3 3 1 3 3 3 3 3 3 1 1 3 3 3 3 1 3 1 3 1 3 3 1 3 3 3 3 3 3 1 1 3 3 3 1 3 3 3 1 3 3 3 1
[148] 3 3 1

Within cluster sum of squares by cluster:
[1] 44.087545  5.163861 47.450194  5.726442  1.974750
 (between_SS / total_SS =  82.5 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"         "ifault"      

Comment déterminer le nombre de clusters ? (1)

Ces méthodes non supervisées, sont sans a priori sur la structure, le nombre de groupe, des données.

rappel : un cluster est composé

Comment déterminer le nombre de clusters ? (2)

Comment déterminer le nombre de clusters ? avec la classification hiérarchique

La coupure de l’arbre à un niveau donné construit une partition. la coupure doit se faire :

plot(iris.scale.hclust.ward, hang = -1, cex = 0.5)

Comment déterminer le nombre de clusters ? avec les kmeans

Comparaison des résultats des deux clustering

29 0 0
20 0 0
1 29 0
0 24 21
0 0 26

Pros et cons des différents algorithmes

Algorithme Pros Cons
Hiérarchique L’arbre reflète la nature imbriquée de tous les sous-clusters Complexité quadratique (mémoire et temps de calcul) \(\rightarrow\) quadruple chaque fois qu’on double le nombre d’individus
Permet une visualisation couplée dendrogramme (groupes) + heatmap (profils individuels)
Choix a posteriori du nombre de clusters
Algorithme Pros Cons
K-means Rapide (linéaire en temps), peut traiter des jeux de données énormes (centaines de milliers de pics ChIP-seq) Positions initiales des centres est aléatoire \(\rightarrow\) résultats changent d’une exécution à l’autre
Distance euclidienne (pas appropriée pour transcriptome par exemple)

Visualisation des données - coloration par espèces

species.colors <- c(setosa = "#BB44DD", virginica = "#AA0044", versicolor = "#4400FF")
plot(mes.iris, col = species.colors[iris$Species], cex = 0.7)

*** # Supplementary materials

Distance avec R : indice de Jaccard

          v.a       v.b
v.b 0.3333333          
v.c 0.0000000 0.3333333

Comparaison de clustering: Rand Index

Mesure de similarité entre deux clustering

à partir du nombre de fois que les classifications sont d’accord

\[R=\frac{m+s}{t}\]

Comparaison de clustering: Adjusted Rand Index

\[ ARI=\frac{RI-Expected RI}{Max RI -Expected RI}\]

Comparaison de clustering: Rand Index

Mesure de similarité entre deux clustering

à partir du nombre de fois que les classifications sont d’accord

\[R=\frac{m+s}{t}\]

Comparaison de clustering: Adjusted Rand Index

\[ ARI=\frac{RI-Expected RI}{Max RI -Expected RI}\]

Comparaison des résultats des deux classifications

clues::adjustedRand(cluster.hclust5, cluster.kmeans3)
     Rand        HA        MA        FM   Jaccard 
0.7848770 0.4637776 0.4730527 0.6167001 0.4299265 
## Print the complete list of libraries + versions used in this session
sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] pheatmap_1.0.12    vegan_2.5-6        lattice_0.20-38    permute_0.9-5      rgl_0.100.30       RColorBrewer_1.1-2 clues_0.6.1        FactoMineR_2.0     kableExtra_1.1.0   knitr_1.28        

loaded via a namespace (and not attached):
 [1] ggrepel_0.8.1           Rcpp_1.0.3              assertthat_0.2.1        zeallot_0.1.0           digest_0.6.24           mime_0.9                R6_2.4.1                backports_1.1.5         evaluate_0.14           highr_0.8               httr_1.4.1              ggplot2_3.2.1          
[13] pillar_1.4.2            rlang_0.4.4             lazyeval_0.2.2          rstudioapi_0.10         miniUI_0.1.1.1          Matrix_1.2-17           rmarkdown_2.1           splines_3.6.1           webshot_0.5.2           readr_1.3.1             stringr_1.4.0           htmlwidgets_1.5.1      
[25] munsell_0.5.0           shiny_1.4.0             compiler_3.6.1          httpuv_1.5.2            xfun_0.12               pkgconfig_2.0.3         mgcv_1.8-31             htmltools_0.4.0         flashClust_1.01-2       tidyselect_0.2.5        tibble_2.1.3            viridisLite_0.3.0      
[37] crayon_1.3.4            dplyr_0.8.3             later_1.0.0             MASS_7.3-51.4           leaps_3.0               grid_3.6.1              nlme_3.1-142            jsonlite_1.6.1          xtable_1.8-4            gtable_0.3.0            lifecycle_0.1.0         magrittr_1.5           
[49] scales_1.1.0            stringi_1.4.6           promises_1.1.0          scatterplot3d_0.3-41    xml2_1.2.2              vctrs_0.2.0             tools_3.6.1             manipulateWidget_0.10.0 glue_1.3.1              purrr_0.3.3             hms_0.5.2               crosstalk_1.0.0        
[61] parallel_3.6.1          fastmap_1.0.1           yaml_2.2.1              colorspace_1.4-1        cluster_2.1.0           rvest_0.3.5            

… par une projection sur une ACP

par(mfrow = c(1,2))
biplot(prcomp(mes.iris), las = 1, cex = 0.7,
       main = "Données non normalisées")
biplot(prcomp(mes.iris, scale = TRUE), las = 1, cex = 0.7,
       main = "Données normalisées")

Cas d’étude : TCGA Breast Invasive Cancer (BIC)

TP : analyse de données d’expression

Contact: